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We determine the optimal scaling of local-update flat-histogram methods with system size by 
using a perfect flat-histogram scheme based on the exact density of states of 2D Ising models. The 
typical tunneling time needed to sample the entire bandwidth does not scale with the number of 
spins A'^ as the minimal A'^^ of an unbiased random walk in energy space. While the scaling is power 
law for the ferromagnetic and fully frustrated Ising model, for the ±J nearest-neighbor spin glass 
the distribution of tunneling times is governed by a fat-tailed Frechet extremal value distribution 
that obeys exponential scaling. We find that the Wang-Landau algorithm shows the same scaling 
as the perfect scheme and is thus optimal. 
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Monte Carlo methods are well-suited for the simulation 
of large many body problems, since the complexity for a 
single Monte Carlo update step scales only polynomially 
and often linearly in the system size, while the config- 
uration space grows exponentially with the system size. 
The performance of a Monte Carlo method is then deter- 
mined by how many update steps are needed to efficiently 
sample the configuration space. For second order phase 
transitions in unfrustrated systems the problem of "crit- 
ical slowing down" - a rapid divergence of the number 
of Monte Carlo steps needed to obtain a subsequent un- 
correlated configuration - was solved more than a decade 
ago by cluster update algorithms . At first order phase 
transitions and in systems with many local minima of the 
free energy such as frustrated magnets or spin glasses, 
there is the similar problem of long tunneling times be- 
tween local minima. With energy barriers AE scaling lin- 
early with the linear system size L, the tunneling times r 
at an inverse temperature f3 = l/kgT scale exponentially 
with the system size, r ~ exp(/3Ai?) oc exp(const x L). 
Several methods were developed to overcome this tun- 
neling problem, such as the multicanonical method 0, 
broad histograms simulated and parallel tempering 
, and Wang-Landau sampling j^] . The common aim of 
all these methods is to broaden the range of energies sam- 
pled within Monte Carlo simulations from the sharply 
peaked distribution of canonical sampling at fixed tem- 
perature in order to ease the tunneling through barriers. 

Ideally, all relevant energy levels are sampled equally 
often during a simulation, thus producing a "flat his- 
togram" in energy space. Some methods approach this 
goal by variations and generalizations of canonical dis- 
tributions 0, 0, while others 0, discard the notion 
of temperature completely and instead are formulated in 
terms of the density of states. With a probability p{E) 
for a single configuration with energy E, the probability 
of sampling an arbitrary configuration with energy E is 



given as Pe = p{E)p{E), where the density of states p{E) 
counts the number of states with energy E. Upon choos- 
ing p{E) cx l/p{E) instead of p{E) oc exp(— one ob- 
tains a constant probability Pe for visiting each energy 
level E, and hence a flat histogram. Wang and Landau [5j 
proposed a simple and elegant flat histogram algorithm 
that iteratively improves approximations to the initially 
unknown density of states p{E). Once p{E) is determined 
with sufficient accuracy, the Monte Carlo algorithm just 
performs a random walk in energy space. Within two 
years of publication this algorithm has been applied to 
a large number ofproblems d, 111 and extended to 
quantum systems |2| . 

In this Letter we investigate the performance of flat 
histogram algorithms in general, and the Wang-Landau 
algorithm in particular, for three systems for which the 
density of states p{E) is known exactly on finite two- 
dimensional (2D) lattices: the Ising ferromagnet as the 
simplest example, the fully frustrated Ising model as a 
prototype for frustrated systems, and the ± J Ising spin 
glass. For each of these models we construct a perfect 
fiat histogram method by simulating a random walk in 
configuration space where we employ the known density 
of states for these models to set p{E) cx 1/ p{E). 

As a measure of performance we use the average tun- 
neling time T to get from a ground state (lowest energy 
configuration) to an anti-ground state (configuration of 
highest energy) , which is the relevant time scale for sam- 
pling the whole phase space [loj . Since the number of 
energy levels in a c?-dimensional system with linear size 
L scales with the number of spins N = L"^, the tunneling 
time for a pure random walk in energy space is 
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This scaling was in fact found for first order phase tran- 
sitions 0, Q. However, none of the systems we study 
exhibit this scaling. While the ferromagnetic and fully- 
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FIG. 1: Scaling of tunneling times r from the ground state 
to the anti-ground state as a function of system size L for a 
perfect flat histogram method that samples using the exact 
density of states. Shown are results for the ferromagnet (FM) 
and fully frustrated (FF) 2D Ising models with both local 
and N-fold way updates; the inset illustrates the frustrated 
couplings. In all cases polynomial scaling r oc 1^"+' is found, 
with = 0.743 ± 0.007, z"' ' 
1.727 ± 0.004, and zjfif^t 
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iv-foid = 0.729 ±0.011, 
1.692 ± 0.004. 



frustrated models exhibit power law scaling, for the spin 
glass the distribution of characteristic tunneling times 
is extremely broad and appears to diverge exponentially 
with system size. For all three models the scaling of the 
performance of the Wang-Landau algorithm is the same 
as that of the perfect flat histogram method. 

We first look at homogeneous systems, with both fer- 
romagnetic and fully frustrated couplings (see inset of 
Fig. nj . Exact densities of states p{E) were calculated 
using the program of Beale fo'^ the ferromagnet and 
the algorithm of Saul and Kardar for the fully frus- 
trated Ising model. In Fig. ^ the measured tunneling 
times are plotted versus system size L for the perfect 
flat histogram method using both local and A''-fold way 
updates. Instead of the expected scaling Eq. we 
find a more rapid increase following r oc L^d+z ^[^]^ 
z™i = 0.743 ± 0.007 for the ferromagnetic (FM) model 
and zll^^ = 1.727 ± 0.004 for the fully frustrated (FF) 
model. Frustration significantly increases the scaling ex- 
ponent but still conserves power law scaling. 

A previous study suggests that TV-fold way up- 
dates 01 speed up Wang-Landau sampling. We find 
identical scaling exponents for A'^-fold way and local up- 
dates within our error bars, see Fig.^ implying that any 
performance improvement remains constant with system 
size. As can be seen from Fig. 12 A^-fold way updates 
reduce the tunneling time by roughly a factor of 2, in- 
dependent of system size. In practice, the reduction of 
tunneling times is offset by the added expense of iV-fold 
way updates. The slight decrease of tunneling ratios in 
the fully frustrated model compared to the ferromagnetic 
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FIG. 2: Scaling of the ratios of tunneling times measured 
using local and A^'-fold way updates. All simulation are done 
with single-spin flip updates weighted by the exact density of 
states. For each system size of the ferromagnet and the fully 
frustrated 2D Ising models 50,000 tunneling events and at 
least 25 tunneling events for the 2D spin glass 1000 randomly 
generated realizations for L = 6, 8 ... 18 and 350 realizations 
for L= 20 are used for averaging. 
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FIG. 3: Normalized distributions of tunneling times (left pan- 
els) and the ratio of the number of flrst excited states to the 
number of ground states p{Ei)/ p{Eo) (right panels). For both 
system sizes, L — 6 and L = 16, 1000 randomly generated 2D 
±J spin glass realizations are sampled. The measured his- 
tograms of tunneling times (using local updates and the ex- 
act density of states) and p[E\) / p{Eo) both follow fat tailed 
Frechet distributions (solid lines). 



model can be explained by the increasing degeneracy of 
ground states in the fully frustrated model, which de- 
creases the ground state lifetime and makes A''-fold way 
updates less effective. 

To determine the performance in more complex en- 
ergy landscapes, we study the 2D ±J Ising spin glass, 
where we find exponential scaling. We measured tunnel- 
ing times of the perfect flat histogram method for 1000 
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FIG. 4: Scaling of the parameters of Frechet distribution of 
the tunneling times (see Eq. |2l of the 2D ± J Ising spin glass 
using perfect flat histogram sampling as a function of system 
size L. Solid and dashed lines show least square fits assuming 
exponential and algebraic (power law) scaling respectively. 
The exponential fits of the data with L > 8 for /i and /3 give 
= 12 and 7.6, respectively, being significantly better than 
algebraic fits with = 130 and 61, respectively. The inset 
shows the scaling of the shape parameter. 



realizations for the system sizes L = 6,8, 10, 12, 14, 16, 18, 
and 350 realizations for L = 20 using the exact density 
of states obtained by the algorithm of Saul and Kardar 
[l^ . For fixed system size the tunneling times are scat- 
tered over several orders of magnitude for the various 
realizations, as shown in Figs. Ola) and 13 b). To analyze 
the underlying distributions we use extremal value theory 
[l4|. The central limit theorem for extremal values [13 
states that the extrema of large samples are distributed 
according to one of only three distributions, depending 
on whether the tails of the original distribution are fat- 
tailed (algebraic), exponential or thin-tailed (decaying 
faster than exponential). This theorem is successfully 
applied in the analysis of tails in diverse fields such as 
hydrology, insurance and finance 0|. Surprisingly, here 
we find that not only the extrema, but all of the measured 
tunneling times [see Figs. |21a) and|3|b)] are distributed 
according to the Frechet extremal value distribution for 
fat-tailed distributions: 



Hi-t,-p{T) = exp 
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with ^ < 0. The parameters of the distribution are deter- 
mined by a maximum likelihood estimator. Fig. 0] shows 
that the location parameter specifying the maximum 
of the distribution and the scale parameter [3 determin- 
ing the width of the distribution scale exponentially with 
linear system size L: 



H oc exp(L/(4.21 ± 0.04)) 
/3 oc exp(L/(3.37±0.05)) 
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FIG. 5: Correlation between tunneling times and ratios of the 
number of first excited states to the number of ground states 
p{Ei)/p{Eo). Shown are data from 1000 randomly generated 
2D ± J spin glass realizations of fixed system size L = 16. 



The shape parameter ^, shown in the inset of Fig. 0] 
determines the power law decay of the fat tails of the 
distribution 
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dr 
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From this asymptotic behavior one can see that the m-th 
moment of a fat tailed Frechet distribution (with ^ < 0) 
is weh defined only if |^| < 1/m. We find that |^| > 1/2 
for L > 10, which implies that the variance (m = 2) does 
not exist and the central limit theorem for mean values 
does not apply. Any direct estimate of the mean tun- 
neling time - as opposed to the most likely time given 
by the Frechet location parameter - then has an infinite 
error. This breadth may explain differences in our con- 
clusion from those in Refs. llTl It also looks plausi- 
ble, although we cannot go to large enough systems, that 
1^1 increases monotonically with system size and could 
become larger than 1, in which case even the mean tun- 
neling time (m = 1) becomes ill-defined. The most likely 
tunneling time, given by the location parameter /i, would 
still remain well defined and finite. 

Since the tunneling time is to a large extent domi- 
nated by the energy landscape at low energies, we study, 
as a qualitative measure for the complexity of the energy 
landscape, the distribution of p{Ei) / p{Eo) where p{Ei) 
is the number of first excited states and p{Eq) the num- 
ber of ground states. Again we find fat tailed Frechet 
distributions, as shown in Figs. |31c) and|21d), suggesting 
that intrinsic properties of the 2D ± J spin glass account 
for the observed distribution of the measured tunneling 
times. In Fig. we show the measured tunneling times 
versus the ratio p{Ei)/p{Eo). A strong correlation over 
five orders of magnitude is found, supporting the argu- 
ment. 
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FIG. 6: Convergence of the tunneling times r during the cal- 
culation of the density of states using Wang Landau sampling 
versus Wang-Landau parameter /. Main panel: 2D Ising fer- 
romagnet model; inset: typical samples for the 2D ±J spin 
glass. For each system size the tunneling times are averaged 
over 50 (FM) and 500 (SG) independent runs, and are given 
in units of the average tunneling time Toxact measured dur- 
ing random walks in the exact density of states using Wang- 
Landau (local update) sampling. 

The question arises why the scaling behaviors of the 
fully frustrated model and the spin glass are differ- 
ent. Both models have an exponentially large number 
of ground states and an extensive ground state entropy, 
but the tunneling time scaling is algebraic for the fully 
frustrated model and exponential for the spin glass. We 
believe that the reason is the difference in complexity of 
the energy landscapes in the two models. While the en- 
ergy landscape above the large number of ground states 
of a frustrated model can be simple, the energy landscape 
of the spin glass is more complex with an extremely large 
number of local minima: The number of first excited 
states p{Ei) that can be reached from the p{Eo) ground 
states by one single spin flip is at most L'^p{Eo). The 
number of local minima that are not connected to the 
ground state by a single spin flip can thus be estimated 
as p{Ei) — L'^p{Eq). For the 2D ferromagnetic model the 
ratio p{Ei)/p{Eq) is exactly L^, and p{Ei) / p{Eq) is of 
order for the fully frustrated model. In contrast, in 
the 2D spin glass the ratio p{Ei) / p{Eo) exceeds by 
several orders of magnitude for some realizations, as can 
be seen from Figs. |21c),|31d) and[Sl 

Finally we compare the tunneling times measured in 
the Wang-Landau algorithm to the perfect flat histogram 
method. The Wang-Landau algorithm approaches the 
exact density of states p{E) by multiplying the current 
estimate at each visited level by a factor / that is re- 
duced towards 1 over time as the algorithm converges. In 
Fig. El we show tunneling times for Wang-Landau sam- 
pling as a function of this correction factor / for the Ising 
ferromagnet (main panel) and for the spin glass (inset). 
Results for the fully frustrated model (not shown) are 



qualitatively similar. In the initial stages of the simula- 
tion (In / ^ 10~^) the tunneling times are shorter than 
for exact sampling, since the random walk is biased — 
it is always driven away from the last region visited (due 
to the increased p{E) there). Eventually the tunneling 
times converge to exactly the same times as for the per- 
fect flat histogram method, indicating convergence of the 
Wang-Landau algorithm. The Wang-Landau algorithm 
is thus optimal in the sense that it performs identically to 
a perfect flat histogram method. Unlike in recent applica- 
tions to continuum systems (19j no convergence problems 
are observed for these lattice models. 

Our benchmarks of a perfect flat histogram method 
provide a lower bound for the tunneling times of other 
" flat histogram" methods such as multicanonical I3| , tem- 
pering0], broad histograms 0| or Wang-Landau sam- 
pling |5(. From our analysis of tunneling times we flnd 
that the Wang Landau algorithm scales identically as 
the perfect flat histogram method and is thus optimal. 
We expect that the other methods will perform simi- 
larly when well-tuned. However, whether one uses local 
or iV-fold way updates, the scaling is not the iV^ scal- 
ing of a random walk in a system with N energy lev- 
els, but slower, namely N'^L^ for both the ferromagnetic 
(z = 0.743 ± 0.007) and the fully frustrated Ising model 
(z = 1.727 ± 0.004). The power law scaling for the frus- 
trated model is very encouraging and demonstrates, for 
the first time, that the Wang-Landau algorithm is well 
suited for frustrated models. A combination with alter- 
native sampling schemes, such as in Ref. jS^j can further 
improve the performance. The observation of a power 
law with an exponent larger than 2 is not trivially ex- 
plained in the context of a random walk in energy space 
and the subject of further investigations. 

The exponential scaling for the ± J spin glass even for 
the perfect method shows a limitation for any flat his- 
togram method. Here the distribution of tunneling times 
follows a fat tailed Frechet extremal value distribution. 
The origin of this extremal character of the 2D ± J Ising 
spin glass remains an interesting open question. Fur- 
ther studies are in progress to investigate this issue as 
well as three-dimensional classical spin and quantum spin 
glasses. 
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